function eigen

%  Solves the eigenvalue problem:
%       [p(x)y']' - r(x)*y = - l^2*q(x)*y   for xL < x < xR
%  where
%      y(xL) = 0  and y(xR) = 0

clear *
%format short e


% set endpoints
	xL=0; 
	xR=1; 

% set number of points along the x-axis
n=10000

% generate the points along the x-axis: x(1)=xL and x(n+2)=xR
x=linspace(xL,xR,n+2);
h=x(2)-x(1);

%%%%%%%%%  symmetric approximation
% calculate the coefficients of finite difference equation
a=zeros(1,n); b=zeros(1,n); c=zeros(1,n);
for i=1:n
	hL=x(i+1)-h/2;
	hR=x(i+1)+h/2;
	pL=p(hL);  pR=p(hR);  qI=q(x(i+1));  rI=r(x(i+1));
	a(i)=(pR+pL+h^2*rI)/(h^2*qI);
	b(i)=-pL/(h^2*qI);
	c(i)=-pR/(h^2*qI);
end;

Am=zeros(n,n);
for i=1:n
	Am(i,i)=a(i);
	if(i<n)
		Am(i,i+1)=c(i);
	end
	if(i>1)
		Am(i,i-1)=b(i);
	end
end

%Am

es = eigs(Am,7,'sm')

function f=p(x)
f=1;

function f=r(x)
f=0;

function f=q(x)
%f=(3*pi/7)^2*(1+x)^4;
f=(1+x)^4;

% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end
